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IX 


1.  INTRODUCTION 


Many  circumstances  exist  for  which  one  may  wish  to  identify,  or  at  least  image,  a 
distant  object  (target),  for  example,  a  vehicle,  when  one  cannot  see  the  object  clearly 
because  of  obscuration,  such  as  fog,  smoke,  or  haze.  A  useful  technique  for  performing 
such  imagery  may  be  laser  imaging  radar,  or  lidar.  Key  advantages  of  lidar  are  (1)  using 
short  pulses,  one  may  detennine  a  range  profde,  (2)  using  receiving  optics,  one  may 
determine  azimuth  and  elevation  resolution.  Because  the  transmitted  electromagnetic 
radiation  is  either  visible  or  infrared,  the  relatively  short  wavelength  provides  a 
diffraction-limited  angular  resolution  much  finer  than  for  microwave  radar,  for  which 
target  rotation  is  usually  required  to  obtain  fine  cross-range  resolution. 

A  particularly  useful  lidar  concept  is  “flash  lidar.”  A  single  laser  pulse  illuminates 
the  entire  target.  Scattered  light  (photons)  from  various  parts  of  the  target  is  received  by 
the  lidar  optics  and  focused  onto  an  array  of  detectors.  Because  the  target  is  usually  quite 
far  from  the  lidar,  the  array  of  detectors  (called  the  “focal  plane  array,”  or  FPA)  is 
typically  in  the  focal  plane.  A  target  image  is  thereby  formed;  each  detector  produces  a 
single  pixel.  Each  detector  is  also  assumed  capable  of  resolving  the  arrival  time  of  a 
photon.  When  this  arrival  time  is  interpreted  as  a  delay  td  relative  to  the  arrival  time  of  a 
photon  from  the  part  of  the  target  nearest  the  lidar,  then  the  depth  of  the  pixel  may  be 
inferred  as  cf/2,  and  a  3D  image,  consisting  of  a  3D  array  of  voxels,  may  be  formed. 

Some  photons  travel  through  the  obscurant  without  scattering;  these  are  termed 
“ballistic”  photons  and  can  produce  an  undistorted  3D  image.  Others  are  scattered  but 
nevertheless  reach  the  FPA;  these  “diffuse”  photons  interfere  with  the  clarity  of  the  target 
image.  For  given  assumptions  about  the  parameters  of  the  lidar  (including  detector  noise 
level),  obscurant,  and  target,  our  model  simulates  a  3D  target  image  and,  for  each  voxel, 
indicates  the  number  of  detector  counts  due  to  ballistic  photons,  diffuse  photons,  and 
noise. 
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2.  FLASH  LID AR 


We  assume  a  circular  receive  aperture  (lens  or  mirror)  of  diameter  D  and  focal 
length  d,  with  an  FPA  consisting  of  nx  x  nv  detectors  (pixels),  with  center-separation 
(pitch)  w«d  .  The  single-pixel  (“instantaneous”)  field-of-view  (IFOV)  is  then  w/d 
(radians)  and  the  overall  field-of-view  (FOV)  is  ~  nx  w/d  x  ny  w/d.  Because  the  target 
and/or  lidar  platform  may  be  moving,  it  is  advantageous  to  transmit  a  short  pulse  and 
receive  the  echo  pulse  on  all  detectors  in  parallel.  Each  detector  records  the  arrival  time 
of  photon(s),  quantizing  this  into  a  series  of  “time  bins.”  This  concept  is  “flash  lidar”  as 
opposed  to  the  more  traditional  scanning  lidar.  The  plane  in  which  the  target  is  focused  is 
the  image  plane.  We  assume  that  the  target  is  relatively  far  from  the  lidar,  and  that  the 
image  plane  is  therefore  the  focal  plane.  We  also  assume  that  the  flash  lidar  is  not  a 
search  device,  but  that  the  target  has  been  already  located  by  another  sensor,  for  example, 
a  microwave  radar  or  a  forward-looking  infrared  (FLIR).  The  lidar  is  a  “soda-straw” 
sensor,  used  for  target  recognition  and  identification.  Figure  1  presents  an  overview  of 
the  flash  lidar  concept. 

Outgoing  Laser  Pulse 
(Beer’s  Law) 

Isotropically  Scattering 
Spherical  Target 
Return  Path  (Monte  Carlo) 

“Absorption  Gated”  Photon 
Aperture 

“Spatially  Gated”  Photon 
Focal  Plane  Array 


“Diffuse”  Photon 


“Ballistic”  Photon 


Figure  1.  Flash  Lidar — Overview 
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The  lidar  transmits  a  pulse  of  energy  E  (millijoules)  and  width  z  (nanoseconds). 
The  beam  is  assumed  to  be  well  collimated,  such  that,  in  the  absence  of  obscurants,  most 
of  the  energy  strikes  the  target,  which  is  at  range  R ,  and  all  of  the  target  is  illuminated  by 
the  beam.  Scattering  from  the  target  is  assumed  to  be  isotropic  (see  Section  4).  Scattered 
energy  is  received  by  the  lidar  optics  and  focused  onto  the  FPA. 

If  obscurant  is  present  (assumed  for  now  to  be  homogeneous  everywhere  in 
space),  then  lidar  energy  moving  a  distance  z  through  the  obscurant  will  decay  according 
to  exp(-z/Z),  where  z  is  the  distance  traveled,  and  L  is  the  scattering  length  and  equals  the 
mean  free  path  of  the  photon  ([1],  p.  259).  This  simple  equation  is  known  as  Beer’s  Law 
([2],  Vol.  3,  p.  4).  The  unscattered  photons  continue  to  travel  in  a  straight  line.  After 
Wang  et  al.  [3]  we  refer  to  these  as  “ballistic”  photons.  Each  time  a  photon  is  scattered, 
the  probability  of  its  not  being  absorbed  is  the  single-scattering  albedo;  the  photons  that 
are  scattered  but  not  absorbed  are  referred  to  as  “diffuse”  photons.  (Wang  et  al.  [3]  also 
refer  to  “snake”  photons  as  those  that  are  scattered,  but  not  enough  to  hit  a  different 
detector  element.  We  prefer  to  avoid  this  term.)  Values  of  y=\!L  for  some  typical 
obscurants  and  wavelengths  are  given  in  [4]. 

When  a  fog  or  other  atmospheric  obscurant  has  a  “visibility”  characterized  by  a 
particular  range,  then,  at  that  range,  the  received  intensity  is  two  percent  of  the  clear-air 
value.  Thus  the  “visibility”  =  3.912  scattering  lengths  [4],  since  e~3'912  =  0.02. 


When  a  pulse  is  transmitted,  the  estimated  number  of  ballistic  photons  received 
by  a  single  detector  (pixel)  is  then 


N 


b 


E  7z(D/2)2  1 

- p - ; - 11  -e 

{he!  A)  H  AnR 2  Npix 


(1) 


where  A  is  the  photon  wavelength,  h  is  Planck’s  constant,  c  is  the  speed  of  light,  he/ A  is 
the  energy  of  a  single  photon,  E  is  the  pulse  energy,  R  is  the  range  to  the  target,  L  is  the 
scattering  length,  D  is  the  aperture  diameter,  p  is  the  target  reflectivity,  Npix  is  the  number 
of  pixels  occupied  by  the  target  image  in  the  FPA.  77  is  the  “overall  efficiency”  and 
includes  such  parameters  as  overlap  of  beam  and  target,  optics  transmission,  FPA  fill 
factor  (ratio  of  area  occupied  by  detectors  to  total  FPA  area),  and  detector  quantum 
efficiency.  We  assume  that  7  =  0.1.  We  assume  that,  by  the  time  the  first  ballistic  photons 
arrive  at  the  FPA,  the  rate  of  arrival  of  received  “pathback”  photons  (those  scattered  from 
the  obscurant  as  the  pulse  travels  to  the  target)  is  small  compared  with  the  rate  of  arrival 
of  ballistic  and  diffuse  photons  scattered  from  the  front  (nearest  point)  of  the  target.  This 
assumption  is  supported  by  experiments  [5]. 
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Before  performing  a  detailed  Monte  Carlo  calculation,  it  is  important  to  compute 
Nh-  If  Nb  is  large,  say  >  1,000,  then  it  is  probably  straightforward  to  obtain  a  clear  image. 
If  Nb  is  less  than  ~1,  then  a  clear  image  may  not  be  possible  using  a  single  laser  pulse, 
though  it  may  be  possible  using  pulse  integration  and  a  very  low-noise  detector.  Figure  2 
shows  Nb  for  some  specific  parameter  values. 


h:- 6.626 10  ^-J-s  c:=3-10^-m-s 

X:=  1.54 10_  6-m  Epulse  :=4010_  ’  J 


D  :=  6- cm  ^overall  :=  0. 1  Nx:=27  Ny:=67  ptgt  :=  0. 1 


Epulse  -  RoverL:  l  2 )  1 

Nph.  .  :=  — - — e  - ptet-ri  overall 

bJ  4k^2  Nx-Ny 


NBph.  . 


Nph. 


.•e 

J 


RoverL  j 


R  over  L 

b™  R  =  0.1  km 
R  =  0.32  km 
KHX  R=  1  km 


R  over  L 

DnD  R  =  0. 1  km 
--o-  R  =  0.32  km 
XXX  R  =  1  km 


Figure  2.  Estimated  Number  of  Total  and  Ballistic  Photons  per  Pixel 
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3.  PHOTON  SCATTERING 


To  consider  the  scattering  of  a  photon,  we  use  two  spherical  coordinate  systems  - 
a  “lab  frame,”  with  the  z-axis  along  the  initial  photon  direction,  and  a  “photon  frame,” 
whose  z  axis  is  along  the  current  direction  of  propagation  (initially,  the  two  frames  are 
co-aligned).  At  each  scattering  event  the  new  direction  and  new  photon  frame  is 
characterized  by  polar  angle  9  and  azimuthal  angle  (f>,  measured  with  respect  to  the 
previous  photon  frame.  The  distribution  of  scattering  directions  is  given  by  the  phase 
function  P(0,<f>)  ([6],  Vol.  2,  Chapter  10).  We  assume  a  simple  model  such  that  (f)  is 
distributed  unifonnly  from  0  to  2n  (photon  polarization  effects  are  neglected),  and  the 
distribution  over  9  is  given  by  0(6)  (the  standard  notation  [6]).  The  nonnalization  is  as 
follows: 

P{9,(f>)d9d(j)  =  Q>(9)sm{9)d9d(j) 

Inn 

J  j  <I>(0)sm(0)d6)d</>  =4/r  .  (2) 

0  0 

n 

jo(9)sin(9)d9  =2 

0 


The  actual  form  of  O  is  given  by  Mie  theory.  Because  it  is  complicated,  a  simple 
approximation  is  typically  used.  Many  such  approximations  exist.  We  chose  the  Henyey- 
Greenstein  (HG)  function  ([6],  Vol.  2,  p.  307)  because  it  is  widely  used  and 
mathematically  tractable: 


0(g,6) 


_ _  ,0<g<  1 

(l  +  g2  -2gcos6>)3  2 


(3) 


The  parameter  g  characterizes  the  degree  of  forward  scattering.  Isotropic  scattering 
occurs  for  g  =  0.  Plots  of  O (g,9)  are  given  in  Figure  3.  The  procedure  for  generating  a 
random  9  according  to  this  distribution  is  given  in  Appendix  A.  Approximate  values  of  g 
for  some  obscurants  are  also  given  in  [7],  based  on  [6],  Vol.  2,  Chapter  10. 
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Scattering  Amplitude  for  Various  g 


®{g,0)=- 


l-g2 


,3/2 


(l  +  ^2-2gcos(6>)) 

g(typical,  1.54  ^im)= 


~0.05  -  haze 

~0.6  --  smoke 

~0.8  -  0.9  -  fog,  cloud,  rain 


Forward  Scattering  Angle:  0 


g  =  0  g  =  0.3  g  =  0.6  g  =  0.9 


Figure  3.  Scattering  Angle  Distributions 
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4.  MONTE  CARLO  MODEL 


The  IDA  Flash  Lidar  Monte  Carlo  Model  works  as  follows.  (An  earlier  version 
was  described  in  [7].)  From  the  lidar  to  the  target,  the  intensity  of  the  beam  is  computed 
simply  by  Beer’s  Law.  We  assume  no  absorption  by  the  obscurant  (often  a  good 
assumption  [6]).  The  total  number  of  photons  that  are  scattered  isotropically  from  the 
target  and  return  to  a  detector  in  the  lidar  FPA  is  then  given  by 


N  = 

total 


(i he  /  A) 


P 


n(D!  2)2  1 


4nR 2  N 


„  -R/L  at 

rj-e  -  N  he 


(4) 


pix 


because  each  photon  scattered  from  the  target  eventually  reaches  the  surface  of  the  sphere 
of  radius  R  centered  at  the  target. 

Each  photon  scattered  from  the  target  is  modeled  separately.  The  photon  leaves 
the  target  in  a  specified  direction  (0  =  0.  0  =  0  in  the  lab  frame).  The  path  length  I 
between  successive  scattering  events  is  selected  randomly  from  an  exponential 
distribution  f(I)  =  yLe~l/L;  this  is  carried  out  by  choosing  u  from  a  uniform  random 
distribution  on  [0,1]  and  computing  /  =  -L  ln(w)  (see  also  [7]).  The  new  scattering 
direction  is  then  chosen  by  randomly  selecting  the  scattering  angles  0  and  </>,  as  described 
above.  The  transformation  from  lab  frame  is  accomplished  by  a  suitable  product  of 
rotation  matrices,  the  exact  form  of  which  is  given  in  Appendix  C.  This  scattering 
procedure  is  iterated  until  the  photon  range  exceeds  R ,  at  which  point  the  model 
computes  the  exit  time,  values  of  0,(f),  and  values  of  0,0.  We  then  invoke  the  assumption 
of  isotropic  scattering  from  the  target:  for  each  photon  leaving  the  target  at  0,0  and 
exiting  the  sphere  at  (0,0),  there  is  an  equally  likely  photon  leaving  the  target  at  (n- 
0,-0  )  and  exiting  the  sphere  at  (0,0).  Therefore  each  modeled  photon  can  be  considered 
as  a  photon  that  strikes  the  lidar  aperture.  This  procedure  greatly  decreases  computation 
time.  Computer  time  is  also  saved  by  tenninating  consideration  of  any  photon  whose 
travel  time  is  greater  than  the  time  corresponding  to  the  last  time  bin  being  recorded  by 
the  detector.  The  model  also  includes  the  option  of  assuming  that  the  obscurant  is 
confined  to  a  layer  between  concentric  spheres  of  radii  Ri  and  Rj  (0  <  R/  <  R2  <  R). 

Many  photons  that  strike  the  aperture  still  do  not  strike  the  FPA,  since  their 
incident  angle  is  too  large;  these  are  “spatially-gated”  photons — see  Figure  1.  For  each 
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photon  striking  the  aperture,  whether  ballistic  or  diffuse,  the  model  computes  the  FPA 
detector  element  (nx,  n  v)  that  it  strikes  (if  any),  and  the  delay  time  t. 

Three-dimensional  plots  of  incident  photons  versus  nx  (width),  ny  (height),  and  t 
(depth)  are  produced.  The  time  spacing  At  corresponds  to  the  FPA  time  resolution;  the 
laser  pulse  width  is  assumed  to  be  much  shorter.  The  3D  voxel  spatial  dimensions  are 
then  wxwx  cAt/2  (assuming  a  FPA  fill  factor  of  1 .0). 

The  user  specifies  the  mean  number  of  noise  counts  (“noise  photons”)  per  voxel, 
which  depends  on  the  type  of  detector  being  used.  The  number  of  noise  counts  in  the 
voxel  is  then  chosen  randomly  according  to  a  Poisson  distribution.  The  model  considers 
both  linear  detectors,  which  record  the  number  of  incident  photons  plus  noise  counts,  and 
Geiger  detectors,  which  record  “1”  if  the  total  of  incident  photons  plus  computed  noise 
counts  in  the  voxel  exceeds  unity.  The  model  accounts  for  the  fact  that,  once  a  Geiger 
detector  has  recorded  a  photon,  its  relatively  long  recovery  time  prevents  it  from 
detecting  another  photon  scattered  from  the  same  laser  pulse.  (We  assume  that  the  pulse 
repetition  frequency  [PRF]  is  low  enough  that  the  detector  recovers  before  the  next  pulse 
is  received.)  The  model  also  allows  for  integration  of  a  number  of  laser  pulses.  After 
pulse  integration,  the  mean  of  the  noise  is  subtracted.  It  is  assumed  that  the  sensor  is  well 
characterized,  the  noise  is  stationary,  and  the  mean  of  the  noise  is  well  known. 

We  expect  the  numbers  of  ballistic  and  diffuse  photons  in  a  voxel  to  be  Poisson- 
distributed  (and  we  specifically  assume  the  noise  to  be  Poisson),  since  they  all  satisfy 
criteria  described  in  [9]:  “The  Poisson  distribution  describes  the  population  of  events  in 
any  interval  of  x  (e.g.,  space  or  time)  whenever  (a)  the  number  of  events  in  any  interval  of 
x  is  independent  of  that  in  any  other  non-overlapping  interval;  (b)  in  any  small  Ax,  the 
probability  of  one  event  is  AAx  and  the  probability  of  two  or  more  vanishes  at  least  as  fast 
as  {Ax)  as  Ax  — »  0;  and  (c)  A  does  not  depend  on  x.  Then  [the  mean]  fi  =  Ax;  E{n)  = 
fa ;  Var(n)  =  // .” 

Table  1  presents  a  summary  of  the  inputs  and  outputs  of  the  Monte  Carlo  model. 


10 


Table  1.  Summary  of  Monte  Carlo  Model  Inputs  and  Outputs 

(B  =  ballistic  photons,  D  =  diffuse  photons,  N  =  noise  photons, 

Nmean  =  average  of  noise  counts  over  FPA) 

Given: 

Pulse  energy 
Wavelength 
Range (R) 

Ri,  R2  =  inner  and  outer  ranges  of  obscurant  layer  (target  is  R  =  0) 

(R2  -  Ri)/(L  (L  =  scattering  length) 

Henyey-Greenstein  factor  =  g  (could  be  generalized),  scattering  Albedo  =  a 
Target  reflectivity  (target  isotropic,  “Lambertian,”  or  “specular”) 

Aperture  diameter 
Focal  length 
FPA  size  (nx,  ny) 

Pitch 

Overall  efficiency  (includes  fill  factor,  optics,  transmission,  quantum 
efficiency) 

Time  gate  (assumed  »  laser  pulse  width) 

Number  of  time  gates 
Number  of  pulses  integrated 
Mean  of  noise  (assume  Poisson) 

Detector  type:  linear  or  Geiger 

Compute: 

B,  B  +  D,  B  +  D  +  N,  B  +  D  +  N-  Nmean 
Photons  received  per  voxel  (x,  Ax;  y  ,  Ay;  t,  At) 

3D  plot,  movie 

3D  FPA  plots,  successive  times 
2D  FPA  plots,  successive  times 

2D  FPA  matrices  (numerical  values),  successive  times 

Photons  vs.  angle  and  time 

Photons  vs.  angle,  all  times;  vs.  time,  all  angles 


Figure  4  shows  model  results  giving  the  photon  distribution  versus  delay  time 
(after  the  arrival  of  the  ballistic  photons)  and  scattering  factor  g,  for  a  point  target;  Figure 
5  illustrates  the  photon  distribution  versus  incident  angle  and  g;  and  Figure  6  shows  the 
photon  distribution  versus  delay  time,  incident  angle,  and  g.  In  Figures  4  and  5,  the  color 
scale  indicates  the  fraction  of  photons  in  the  bin.  Each  row  represents  a  histogram  with 
total  value  1.0;  since  there  are  100  rows,  the  total  of  the  color-scale  value  over  the  figure 
is  100.  It  may  be  seen  that,  for  the  parameters  assumed  (R  =  100  m,  R/L  =  4),  for  g  <  0.9, 
the  typical  delay  is  greater  than  1  nsec  and  the  typical  incident  angle  (angle  between 
photon  direction  and  nonnal  to  aperture)  is  greater  than  5  deg.  In  Figure  6,  the  incident 
angle  is  translated  into  “pixels  away,”  that  is,  the  distance  between  the  center  of  the  FPA 
and  the  detector  recording  the  photon,  measured  in  pixel  diameters.  It  may  be  seen  that 
for  g  <  0.9,  most  photons  arrive  more  than  100  pixels  from  the  FPA  center  and  with  a 
delay  time  greater  than  1  nsec.  For  g  =  0.99,  the  opposite  is  true. 
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If  the  beam  divergence  were  small,  then  many  outward-bound  photons,  even 
though  they  had  been  scattered,  could  still  strike  the  target,  thus  increasing  the  estimated 
value  of  Nb  given  by  (1).  Model  results  illustrating  this  effect  are  given  in  Figure  7,  under 
the  assumptions  used  to  compute  the  results  (R  =  100  m,  R/L  =  4,  g  =  0.99;  see  Section 
6).  Our  target’s  70  cm  diameter  subtends  an  angle,  as  seen  from  the  lidar,  of  ±3.5  mrad 
(0.35  m)/(100  m).  Figure  7  shows  that  the  fraction  of  outgoing  beam  that  scatters  but  still 
strikes  the  target  is  ~  5%  of  the  outward-bound  scattered  photons.  About  2%  (=  e~R/L)  of 
the  outward-bound  photons  are  ballistic.  A  correction  to  the  results  in  Section  6  could  be 
made  for  this  effect;  we  have  not  done  so  here. 

Incident  Photon  Distribution:  Delay  Time  after  Ballistic  Return  vs.  g 
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R  =  100  m,  R/L  =  4.  For  g  <_0.9,  Typical  Delay  >_1  nsec. 

Figure  4.  Photon  Distribution  versus  Delay  Time  and  Scattering  Factor 


12 


Delay  Time  after  Ballistic  Return  Delay  Time  after  Ballistic  Return 


Incident  Photon  Distribution:  Incident  Angle  vs.  g 
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R/L  =  4.  For  g  <_0.9,  Typical  Incident  Angle  >_5  deg. 

Figure  5.  Photon  Distribution  versus  Incident  Angle  and  Scattering  Factor 


Distribution  of  Incident  Photons  on  Focal  Plane  Array:  g  *  0 


10  10  10  10  10 
Pixels  Away  from  Target  Pixel 


Distribution  of  Incident  Photons  on  Focal  Plane  Array:  g  =  0.9 


Pixels  Away  from  Target  Pixel 


Distribution  of  Incident  Photons  on  Focal  Plane  Array:  g  =  0.S 


Distribution  of  Incident  Photons  on  Focal  Plane  Array:  g  a  0.99 
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Pixels  Away  from  Target  Pixel 


Figure  6.  Photon  Distribution  versus  Delay  Time,  Pixel  Location,  and  Scattering  Factor 
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Angular  Dispersion  due  to  Obscurant:  g  =  0.99,  Range  =  100m,  R/L  =  4 
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Figure  7.  Photon  Distribution  versus  Beam  Angle  (left  scale  corresponds  to 
jagged  line;  right  scale  corresponds  to  smooth  line) 
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5.  OTHER  CONSIDERATIONS 


5.1  FORWARD  SCATTERING 

The  HG  function  underestimates  the  number  of  scattered  photons  very  near  the 
forward  direction.  More  exact  phase  functions  should  be  used  for  detailed  estimates.  For 
this  reason  we  assume  g  =  0.99  for  the  highly  forward  scattering  considered  herein,  rather 
than  g  ~  0.9,  as  indicated  in  Figure  3. 

The  effect  is  related  to  the  forward-scattering  phenomenon  encountered  in 
classical  electromagnetic  (EM)  scattering  theory.  For  near-isotropic  scattering  from  a 
sphere  of  radius  a,  Jackson  [10],  p.  451,  states: 

The  scattering  in  the  forward  direction  is  a  typical  diffraction  pattern  with 
a  central  maximum  and  smaller  secondary  maxima  [Bessel  function 
pattern],  while  at  larger  angles  it  is  isotropic...  The  total  scattering  cross 
section  is  obtained  by  integrating  over  all  angles.  Neglecting  the 
interference  terms,  we  find... that  the  shadow  diffraction  peak  [width  ~ 

AJ2  m\  gives  a  contribution  of  m  ,  and  so  does  the  isotropic  part.  The  total 

2 

scattering  cross  section  is  thus  2  na  . 

5.2  THERMAL  PHOTONS 

We  estimate  the  number  of  thennal  photons  that  are  received  at  a  detector 
(individual  pixel)  during  a  time  gate.  Because  the  optics  are  presumably  at  approximately 
room  temperature,  we  make  the  approximation  that  the  detector  (area  =  w  ,  neglecting  the 
fill  factor)  is  at  the  center  of  a  hemisphere  of  radius  R  with  an  inner  surface  at 
temperature  300  K  and  emissivity  of  unity.  The  thermal  emission  from  the  hemisphere  is 
given  by  the  well-known  Planck  function 

sr-  ,  (5) 

where  h,  k,  and  c  are  Planck’s  constant,  Boltzmann’s  constant,  and  the  speed  of  light 
respectively,  /  =  c/A=  194  THz,  T  =  300  K,  and  e  =  emissivity  =  1.  The  number  of 
thennal  photons  per  pixel  expected  in  a  time  bin  is  then  independent  of  R : 
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(6) 


Nppixth  =  &  71 . 2L_ .  f  d(j)  [  R'  smOdd -cosO 

(t)  o  o 

_  7iB(f)  -  A f  -St-  rjw 2 

"  M 

where  Af=  bandwidth,  assumed  to  be  630  GHz  (corresponding  to  zU  =  5  nm),  St  =  time 
bin  =  1  ns.  For  the  assumed  lidar,  Nppixth  =  1.3  x  10  8  photons,  entirely  negligible. 


5.3  SOLAR  PHOTONS 


The  intensity  of  sunlight  at  Earth  at  1.54  pm  is  Fsoiar  ~  270  W  m  2  pm  1  ([1 1],  p. 
172).  We  conservatively  neglect  atmospheric  absorption  (it  is  negligible  in  clear  air  — 
[11],  p.  122)  or  any  obscurant.  We  consider  direct  sunlight  falling  on  the  target  and 
undergoing  Lambertian  scattering  into  2;rsteradians.  We  assume  that  the  scattering  plane 
is  the  xy  plane,  and  that  the  sun  and  the  lidar  are  both  near  the  z  axis  at  small  values  of  6. 
In  this  case  the  fraction  of  emitted  energy  that  strikes  the  aperture  is  ~QrcJn,  where  Orcv 
is  the  solid  angle  subtended  by  the  lidar  aperture,  as  seen  from  the  target.  The  estimated 
number  of  solar  photons  striking  an  individual  detector  during  a  time  gate  is 


N, 


solar 


tf^iR-IFOV)2 -AA-St-O^p-Ti 

(t)  n 


(7) 


As  with  (6)  for  thennal  photons,  the  expression  for  Nsoiar  is  also  independent  of  R. 

For  the  parameters  assumed  in  Section  6,  Nsoiar  equals  0.05  photons.  For  the  linear 
detector  we  may  neglect  this,  especially  since  direct  sunlight  is  unlikely  when  obscurants 
are  present.  For  the  Geiger  detector,  a  continuous  illumination  by  solar  photons  can 
reduce  the  effective  number  of  photons  recorded,  especially  for  scatterers  located 
relatively  deep  in  the  target  compared  with  the  front  scatterers. 
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6.  RESULTS 


Monte  Carlo  results  were  computed  for  a  simulated  target  consisting  of  1 1  points 
in  3D  space  arranged  approximately  in  the  shape  of  a  small  “car,”  as  shown  in  Figure  8. 
The  points  may  be  thought  of,  schematically,  as  two  front  wheels,  two  headlights,  one 
hood  ornament,  two  top  comers  of  the  windshield,  two  top  comers  at  the  rear,  and  two 
rear  wheels.  The  coordinate  system  used  is 

z  =  longitudinal  axis  of  target  (“car”)  =  axis  of  lidar  =  line  between  lidar  and 
target) 

x  =  horizontal  axis,  that  is,  the  left  front  and  the  left  rear  “wheels”  have  the  same  x 
value 

y  =  vertical  axis,  that  is,  both  “headlights”  have  the  same  y  value 

To  prevent  one  point  from  shadowing  another,  the  11  locations  are  chosen  such  that  no 
two  points  have  exactly  the  same  values  of  nx  ,nv.  The  time  bin  At  is  taken  to  be  0.5  nsec; 
the  pixel  depth  =  Az  =  cAt! 2  =  7.5  cm  =  Ax  =  Ay. 


Figure  8.  The  Target  Model:  “Car”  Consisting  of  11  Points 
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Figures  9-12  are  presented  as  a  computer-printed  summary  of  parameters  used 
for  the  model  run,  and  resulting  “quad-charts”  showing: 

•  top-left:  ballistic  photons 

•  top-right:  ballistic  +  diffuse  photons  (B  +  D) 

•  bottom-left:  ballistic  +  diffuse  +  noise  photons  (B  +  D  +  N) 

•  bottom-right:  ballistic  +  diffuse  +  noise  photons,  minus  noise  mean  (B  +  D  + 
N  -  Nmean) 

The  detailed  assumptions  are  provided  in  the  computer  output  at  the  top  of 
Figures  9-12.  In  summary,  we  assumed 

•  R  =  100  m 

•  an  optical  depth  =  R/L  =  4  (slightly  greater  than  the  “visibility”  condition  of 
(3.912)1). 

•  Nmean  =100  for  the  linear  detector  (see  [12])  and  5*1  O'7  for  the  Geiger 
detector  [13,  p.  340] 

•  E  =  0.04  J  for  the  linear  detector,  0.0004  J  for  the  Geiger  detector 

•  a  PRF  of  100  Hz  for  the  linear-detector  lidar  and  10,000  Hz  for  the  Geiger- 
detector  lidar;  both  have  average  radiated  power  =  4  watts 

•  integration  times  of  10  and  100  msec. 

•  Geiger  recovery  time  >  10  nsec  and  <  100  psec  ([14],  p.  353) 

Results  are  summarized  as  follows.  Color  scale  represents  number  of  photons 
recorded,  and  varies  with  each  plot. 

•  Figure  9:  Linear  Detector,  10  msec  integration  time  (1  pulse): 

•  Figure  10:  Linear  Detector,  100  msec  integration  (10  pulses) 

•  Figure  11:  Geiger  Detector,  10  msec  integration  time  (100  pulses):  Note  that 

deep  scatterers  are  often  not  as  bright  as  front  scatterers,  since  sometimes  a 
diffuse  photon  from  a  front  scatterer  will  trigger  the  Geiger  detector,  thereby 
precluding  detection  of  a  subsequent  photon  from  a  deeper  scatterer  on  the 
same  pulse. 

•  Figure  12:  Geiger  Detector,  100  msec  integration  time  (1,000  pulses) 

For  the  parameters  assumed,  the  results  indicate  the  following.  For  the  linear 
detector,  integration  over  10  msec  (1  pulse)  is  insufficient  to  produce  a  meaningful  3D 
image,  whereas  integration  over  100  msec  (10  pulses)  does  produce  a  meaningful  image. 
For  the  Geiger  detector,  integration  over  10  msec  (100  pulses)  is  suffucient  to  produce  a 
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meaningful  image,  and  integration  over  100  msec  (1,000  pulses)  produces  a  very  clear 
image. 

If  image  quality  metrics  (IQMs)  are  developed  for  flash  lidar  images  of  targets  in 
obscurants,  we  suggest  that  since  the  effects  of  multiple  target  points  are  important 
(diffuse  photons  from  one  target  point  affect  the  appearance  of  other  target  points),  a 
standard  3D  multipoint  target  be  used.  This  is  analogous  to  the  standard  U.S.  Air  Force 
bar  chart  for  evaluating  electro-optical/infrared  imagery.  Key  inputs  to  IQMs  for  various 
lidars  could  include  (1)  the  full- widths  at  half  maximum  (FWHM)  in  x,  y,  and  z  of  the 
single-point  3D  point-spread  function  (PSF),  caused  by  diffuse  scattering;  (2)  the  signal- 
to-interference-plus-noise  ratio  (SINR);  and  (3)  the  gray-scale  fidelity. 
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Flash  Lidar,  Monte  Carlo  Simulation  Parameters:  Linear  Detectors 

Environmental  Parameters: 

Range (R)  =  100  m 
Scattering  Length  (L)  =  25  in 
R/L  =  4 

Target  Reflectivity  (p)  =  0.1 
2-way  Delay  Time  =  0.66713  ps 
Anisotropy  Factor  (g)  =  0.99 
Scattering  Albedo  (a)  =  1 

Lidar  Parameters: 

Pulse  Energy  =  0.04  J 
Wavelength  =  1.54  jim 
Photons/Pulse  =  3.100966e+017 
Overall  Efficiency  (i))  =  0.1 
Number  of  Pulses  Integrated:  1 
Dwell  Time:  10  ms 

Pulse  Repetition  Frequency  (PRF):  100  Hz 

Detector  Parameters: 

Aperture  Diameter  (D)  =  6  cm 
Focal  Length  =  6.67  cm 
f  number  -  1.1117 
Pitch  =  50  pm 
IFOV  =  749.6252  prad 
Target  x-dim  =  9  pixels 
Target  y-dim  =  9  pixels 
Outgoing  Beam  x-dim  =  32  pixels 
Outgoing  Beam  y-dim  =  32  pixels 

Voxel  Resolution: 

Width  Resolution  =  7.4963  cm 
Height  Resolution  =  7.4963  cm 
Depth  Resolution  =  7.4948  cm  =  0.5  ns 
diffraction  limit  =  0.31313  cm 

6-19-2003  16:51:8.742  Runtime:  0  day(s);  0  hour(s);  0  minute(s);  46.109  secon 


Results  for  Single-Point-Target: 

Photons:  Incident  on  Aperture  =  1248 
Photons:  Expected  Ballistic  =  22.8572 
Photons  Collected:  Ballistic  =  26 
Photons  Collected:  Ballistic  ♦  Diffuse  =  91 
Photons  Collected:  Ballistic  +  Diffuse  +  Noise  =  73C 

Results  for  11  Point  "Car"  Target: 

Photons:  Incident  on  Aperture  =  13728 
Photons:  Expected  Ballistic  =  251.4296 
Photons  Collected:  Ballistic  =  274 
Photons  Collected:  Ballistic  +  Diffuse  =  965 
Photons  Collected:  Ballistic  +  Diffuse  +  Noise  =  736 


Figure  9.  Parameters  and  Results  for  Linear  Detector,  10  msec  Integration  Time  (1  Pulse) 
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Flash  Lidar,  Monte  Carlo  Simulation  Parameters:  Linear  Detectors 

Environmental  Parameters: 

Range  (R)  =  100  m 
Scattering  Length  (L)  =  25  m 

Target  Reflectivity  (p)  =  0.1 
2-way  Delay  Time  =  0.66713  ps 
Anisotropy  Factor  (gl  =  0.99 
Scattering  Albedo  (a)  =  1 

Lidar  Parameters: 

Pulse  Energy  =  0.04  J 
Wavelength  =  1.54  pm 
Photons/Pulse  =  3.100966e+017 
Overall  Efficiency  (g)  =  0.1 
Number  of  Pulses  Integrated:  10 
Dwell  Time:  100  ms 

Pulse  Repetition  Frequency  (PRF):  100  Hz 

Detector  Parameters: 

Aperture  Diameter  (D)  =  6  cm 
Focal  Length  =  6.67  cm 
f  number  =  1.1117 
Pitch  =  50  pm 
IFOV  =  749.6252  prad 
Target  x-dim  =  9  pixels 
Target  y-dim  =  9  pixels 
Outgoing  Beam  x-dim  =  32  pixels 
Outgoing  Beam  y-dim  =  32  pixels 

Voxel  Resolution: 

Width  Resolution  =  7.4963  cm 
Height  Resolution  =  7.4963  cm 
Depth  Resolution  =  7.4948  cm  =  0.5  ns 
diffraction  limit  =  0.31313  cm 

6-20-2003  10:22:38.109  Runtime:  0  day(s);  0  hour(s);  3  minute(s);  4.859  second 


Results  for  Single-Point-Target: 

Photons:  Incident  on  Aperture  =  12480 
Photons:  Expected  Ballistic  =  228.5724 
Photons  Collected:  Ballistic  =  223 
Photons  Collected:  Ballistic  +  Diffuse  =  903 
Photons  Collected:  Ballistic  +  Diffuse  +  Noise  =  72£ 

Results  for  11  Point  "Car"  Target: 

Photons:  Incident  on  Aperture  =  137280 
Photons:  Expected  Ballistic  =  2514.2962 
Photons  Collected:  Ballistic  =  2497 
Photons  Collected:  Ballistic  +  Diffuse  =  9201 
Photons  Collected:  Ballistic  +  Diffuse  +  Noise  =  73£ 


Figure  10.  Parameters  and  Results  for  Linear  Detector,  100  msec 
Integration  Time  (10  Pulses) 
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Flash  Lidar,  Monte  Carlo  Simulation  Parameters:  Geiger  Detectors 


Environmental  Parameters: 

Range (R)  =  100  m 
Scattering  Length  (L)  =  25  m 
R/L  =  4 

Target  Reflectivity  (p)  =  0.1 
2-way  Delay  Time  =  0.66713  ps 
Anisotropy  Factor  (g)  =  0.99 
Scattering  Albedo  (a)  =  1 

Lidar  Parameters: 

Pulse  Energy  =  0.0004  J 
Wavelength  =  1.064  pm 
Photons/Pulse  =  2.14248578541128e+015 
Overall  Efficiency  (i|)  =  0.1 
Number  of  Pulses  integrated:  100 
Dwell  Time:  10  ms 

Pulse  Repetition  Frequency  (PRF):  10000  Hz 

Detector  Parameters: 

Aperture  Diameter  (D)  =  6  cm 
Focal  Length  =  13.34  cm 
f  number  =  2.2233 
Pitch  =  100  pm 
IFOV  =  749.6252  prad 
Target  x-dim  =  9  pixels 
Target  y-dim  =  9  pixels 
Outgoing  Beam  x-dim  =  32  pixels 
Outgoing  Beam  y-dim  =  32  pixels 

Voxel  Resolution: 

Width  Resolution  =  7.4963  cm 
Height  Resolution  =  7.4963  cm 
Depth  Resolution  =  7.4948  cm  =  0.5  ns 
diffraction  limit  =  0.21635  cm 

6-20-2003  11:5:38.109 


Results  for  Single-Point-Target: 

Photons:  Incident  on  Aperture  =  900 
Photons:  Expected  Ballistic  =  15.7923 
Photons  Collected:  Ballistic  =  23 
Photons  Collected:  Ballistic  +  Diffuse  =  85 
Photons  Collected:  Ballistic  +  Diffuse  +  Noise  =  85 

Results  for  11  Point  "Car"  Target: 

Photons:  Incident  on  Aperture  =  9900 
Photons:  Expected  Ballistic  =  173.715 
Photons  Collected:  Ballistic  =  186 
Photons  Collected:  Ballistic  +  Diffuse  =  613 
Photons  Collected:  Ballistic  +  Diffuse  +  Noise  =  61? 


Runtime:  0  day(s);  0  hour(s);  0  minute(s);  29.969  secon 


Figure  11.  Parameters  and  Results  for  Geiger  Detector,  10  msec 
Integration  Time  (100  Pulses) 
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Flash  Lidar,  Monte  Carlo  Simulation  Parameters:  Geiger  Detectors 


Environmental  Parameters: 

Range (R)  =  100  m 
Scattering  Length  (L)  =  25  m 
R/L  =  4 

Target  Reflectivity  (p)  =  0.1 
2-way  Delay  Time  =  0.66713  (is 
Anisotropy  Factor  (g)  =  0.99 
Scattering  Albedo  (a)  =  1 

Lidar  Parameters: 

Pulse  Energy  =  0.0004  J 
Wavelength  =  1.064  pm 
Photons/Pulse  =  2.14248578541128e+015 
Overall  Efficiency  (ii)  =  0.1 
Number  of  Pulses  Integrated:  1000 
Dwell  Time:  100  ms 

Pulse  Repetition  Frequency  (PRF):  10000  Hz 

Detector  Parameters: 

Aperture  Diameter  (D)  =  6  cm 
Focal  Length  =  13.34  cm 
f  number  =  2.2233 
Pitch  =  100  urn 
IFOV  =  749.6252  prad 
Target  x-dim  =  9  pixels 
Target  y-dim  =  9  pixels 
Outgoing  Beam  x-dim  =  32  pixels 
Outgoing  Beam  y-dim  =  32  pixels 

Voxel  Resolution: 

Width  Resolution  =  7.4963  cm 
Height  Resolution  =  7.4963  cm 
Depth  Resolution  =  7.4948  cm  =  0.5  ns 
diffraction  limit  =  0.21635  cm 

6-20-2003  11:15:36.734 


Results  for  Single-Point-Target: 

Photons:  Incident  on  Aperture  =  9000 
Photons:  Expected  Ballistic  =  157.9227 
Photons  Collected:  Ballistic  =  164 
Photons  Collected:  Ballistic  +  Diffuse  =  670 
Photons  Collected:  Ballistic  +  Diffuse  +  Noise  =  67C 

Results  for  11  Point  "Car"  Target: 

Photons:  Incident  on  Aperture  =  99000 
Photons:  Expected  Ballistic  =  1737.1501 
Photons  Collected:  Ballistic  =  1836 
Photons  Collected:  Ballistic  *  Diffuse  =  6207 
Photons  Collected:  Ballistic  +  Diffuse  +  Noise  =  62C 


Runtime:  0  day(s);  0  hour(s);  4  minute(s);  20.016  secon 


Figure  12.  Parameters  and  Results  for  Geiger  Detector,  100  msec 
Integration  Time  (1,000  Pulses) 
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7.  POST-PROCESSING  OF  GEIGER  RESULTS 


The  one  count  maximum  per  pulse  per  pixel  in  a  Geiger  avalanche  photodiode 
(APD)  imaging  array  results  in  two  observed  phenomena  when  multiple  pulses  are 
integrated:  (1)  the  intensity  distribution  of  the  incident  light  does  not  linearly  map  to  the 
output  of  the  detector,  (2)  in  the  3D  image,  the  intensity  of  “later”  voxels  (greater  delay 
time)  can  be  suppressed  by  returns  in  “earlier”  voxels  with  the  same  pixel  coordinates. 


We  may  compare  the  expected  number  of  received  counts  from  linear  and  Geiger 
detectors  as  follows.  The  number  of  ballistic  photons  received  per  pulse  is  Poisson 
distributed  with  mean  pballistic: 


i  ballistic 

p(  AT  \  _  r1 ballistic  -  ^ballistic 

1  \ 1 V  ballistic  )  ^ 


N, 


ballistic  * 


(8) 


For  the  linear  detector 


Me, 


N pulses  M ballistic’  a  Me, 


(9) 


The  output  initiated  by  one  photon  incident  on  a  Geiger  detector  element  is 
assumed  to  be  indistinguishable  from  that  produced  by  several  incident  photons. 
Therefore,  the  output  from  a  Geiger  detector  element  for  one  pulse  is  reduced  to  a 
Bernoulli  trial  where  “success”  (a  count  in  this  case)  can  be  defined  as  1  or  more  ballistic 
photons  received  and  we  denote  the  probability  of  that  event  by  p.  The  distribution  is 
binomial  with  mean  np  and  variance  np(  1  -  p),  where  n  =  Npuises: 

P  =  p(N banish)  =  \-e^ 

Mcou„,s=n PuisJl-  e'tec)  •  0°) 

fY  ^  —  i  #  sy— ^ballistic 

U  —  r1 counts  e 


Figure  13  illustrates  expected  counts  for  ideal  linear  and  Geiger  detectors.  The  output 
from  the  Geiger  detector  increasingly  diverges  from  that  of  the  linear  detector  as  p ballistic 
increases. 
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Output  (counts)  from  Linear  and  Geiger  Detectors 


Average  Number  of  Incident  Photons  per  Pulse 


Figure  13.  Output  (Counts)  from  Linear  and  Geiger  Detectors 

In  order  to  study  the  effect  of  the  difference  in  output  between  the  two  detector 
classes,  consider  an  image  with  two  pixels:  Br  and  Dm.  Pixel  Br  is  brighter  than  pixel  Dm 
when  the  image  is  obtained  in  linear  mode,  by  a  ratio  known  as  the  contrast  ratio 
CR(Br,Dm).  For  the  linear  detector, 

OR,  (Br,  Dm)  =  /icoui"*(Br\  =  ^ons{Br)  (H) 

Mounts  iDm)  Mphotons  (Dm) 


For  a  Geiger  detector,  the  measured  contrast  ratio  between  the  same  two  regions  is: 


CRg  ( Br,  Dm)  = 


Me, 


t  D  \  ( 1  —  pMPho,om^Br^  \ 

( Br )  l1  ^  ) 


M counts  {^Dni) 


(12) 


Figure  14  shows  the  relationship  between  various  values  of  CR(,  and  jUpiwtons(Br). 
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Contrast  Ratio  (CR)  Degradation  As  a  Result  of  Imaging  in  Geiger-Mode 


io'2  10'1  10°  101  io2 

Number  of  Incident  Photons  per  Pulse  in  Brighter  Region  (Br) 


Figure  14.  Geiger  Detector:  Contrast  Ratio  of  Image 


To  help  remedy  the  degradation  in  CR  resulting  from  the  detector  saturation 
described  by  (12),  a  Geiger  image  may  be  remapped  to  a  scale  more  proportional  to  the 
number  of  incident  photons.  Solving  (10)  for  //, ballistic  and  multiplying  by  NpuiSes,  we  have 

f  \ 


^ remapped  ^ pulses  M ballistic  ^ pulses 


Me, 


N 

y  pulses  y 


(13) 


We  presumably  have  only  one  image  and  therefore  do  not  know  //counts  for  each  pixel;  we 
thus  replace  //COunts  with  our  best  guess,  /Vcounls.  Figure  15  illustrates  this  Geiger  remapping 
technique. 
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Notional  Geiger-Mode  Gray-Scale  Output 


1 


r 


N  =  -N„.  In 


pulses 


N 


counts 


1 


N 

\  pulses  J 


Figure  15.  Geiger  Remapping  Technique 

Applying  (13)  to  a  3D  image  requires  accounting  for  the  recovery  time  needed  by 
the  detector  elements  after  being  triggered.  For  simplicity,  we  assume  the  recovery  time 
Tr  corresponds  to  a  greater  depth  than  the  maximum  depth  to  which  the  target  is  being 
imaged.  This  latter  depth  corresponds  to  some  number  of  time-bins  n,.  In  this  case,  as 
soon  as  an  element  is  triggered  for  a  voxel  located  at  Xuyptk,  the  voxels  “behind”  that  one 
(v;,y/,4+/. .  .xhyj,tnt)  will  be  unresponsive  to  any  incident  photons  arriving  from  the  target 
resulting  from  the  current  laser  pulse.  Therefore,  when  remapping  the  counts 
accumulated  in  a  given  voxel,  a  better  value  of  Npuises  to  use  in  (13),  which  will  be  called 
the  effective  number  of  pulses  integrated  or  Npuises*,  is  simply: 

k- 1 

N*PuiseS(xi’yjdk)=Npulses-YJNcounts{xi’yjdi)  •  (14) 

i=i 

N counts  (xu  Vj,  ti)  is  either  0  or  1.  The  correction  is  not  exact,  since  a  “1”  may  sometimes 
indicate  more  than  one  count.  However,  the  remapping  technique  is  quite  useful.  Figure 
16  depicts  three-dimensional  images  before  and  after  remapping  and  shows  that  the 
results  of  Geiger  remapping  can  result  in  an  image  that  is  much  easier  for  a  human  to 
interpret  than  the  corresponding  raw  (non-remapped)  image. 
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Figure  16.  Three-Dimensional  Geiger  Images:  Raw  and  Remapped 

Finally,  if  a  detector  element  at  coordinate  x„y,  in  the  imaging  array  is 
continuously  illuminated  (e.g.,  by  solar  radiation)  at  a  level  of  p photons  per  time  bin, 
whether  or  not  that  element  is  available  to  be  triggered  in  voxel  xbyp  4  depends  on 
whether  or  not  it  was  triggered  in  any  of  the  preceding  voxels:  x;-,y/,4,  ...,xbyj,tk-i-  This  is 
equivalent,  once  again,  to  a  Bernoulli  trial  where  “success”  can  be  defined  as 
“untriggered.”  The  probability  that  a  detector  element  is  untriggered  at  4  is: 

Puntriggered  (*P  T  j  ,h)  =  P  (N photom  N photons  (Xi  >  X j  P k-l)  =  0)  »  (15) 

and  P(Nphotons(Xi,yj,ti )  =  0)  is  the  “q  ”  ( I  -p)  to  the  p  given  in  (12).  Therefore, 

Puntriggered  (*/ .  Xj .  h  )  =  V  '  =  e  .  (16) 

The  average  effective  number  of  pulses  integrated  for  voxel  xbyp  4  is 

N"P*„.(xby,.h)=Nr^el'»-{l-')  .  (17) 

N pUiSes*  under  various  levels  of  continuous  illumination  is  shown  in  Figure  17.  The 
average  number  of  counts  expected  from  a  portion  of  the  target  located  in  voxel  xhyptk 
then  decreases  along  with  Npuises*- 
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Suppresion  due  to  Continuous  Illumination 


Figure  17.  Geiger  Detector:  Eclipsing  Effect  in  Late  Time  Bins 
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8.  COMPARISON  WITH  DIFFUSION  THEORY 


In  this  section  we  discuss  a  comparison  of  the  Monte  Carlo  results  with  results 
from  diffusion  theory.  The  following  summary  is  for  all  scattering  angles,  for  RJL  »  1, 
which  we  compare  to  our  Monte  Carlo  results. 


8.1  DIFFUSION  EQUATION— PARTICLE  FLOW 

The  simple  diffusion  equation  for  particle  flow  is  obtained  as  follows  ([15],  p.  4). 

2  —1 

The  particle  flow  (or  probability  current  density)  q  (m~  s  )  is  given  by 

q=~DVp  ,  (18) 

where  D  is  the  diffusion  coefficient  (nr  s  )  and  p  is  the  particle  density  (m  ).  From  the 
equation  of  continuity 


Y7  dP 

V -q  = - 

dt 


(19) 


we  have 


W =-^ 
d  dt 


(20) 


which  is  the  diffusion  equation. 


8.2  GREEN  FUNCTION  SOLUTION 

For  this  theoretical  discussion  we  take  t  =  0  as  the  time  the  pulse  leaves  the 
(point)  target.  For  an  infinite  uniform  medium  and  an  impulse  at  r  =  0,  t  =  0, 
corresponding  to  the  laser  pulse  just  leaving  the  target,  the  mathematical  conditions  are 

p0(r,°)  =  ^(r),  p0(cc,t)  =  0  ,  (21) 


and  the  solution  ([15],  p.  29;  [16],  p.  868)  is 


Po(Pt) 


1 

(AnDtf1 


-r2 /4DI 

e 


(22) 


The  solution  is  simply  a  spherically  symmetric  Gaussian  with  a  variance  that  grows 
linearly  with  time;  (22)  is  also  known  in  some  communities  as  the  Smirnov  density  ([17], 
Vol.  1,  p.  206).  Normalization  is 
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(23) 


4ft  oo 

J  JQj dr-p0(r,t )  =1 

o  o 

For  isotropic  scattering,  D  =  cL/3,  where  L  is  the  scattering  length  =  mean  free  path, 
while  for  non-isotropic  scattering,  D  =  cL/( 3(l-<cos#>))  [18],  This  relation  for  the 
diffusion  constant  is  derived  in  Appendix  C.  For  the  Henyey-Greenstein  approximation, 
it  is  easily  shown  that  <cos#>  =  g,  where  <  >  indicates  an  ensemble  average  over  the 
distribution.  For  long  times,  the  density  decays  as  t  3/2. 


8.3  TOTAL  PARTICLE  EXIT  RATE  PAST  r  =  R 

The  expected  number  of  particles  in  a  sphere  of  radius  R  at  time  t  is  given  by 


A 

P(R,t)  =^dr-4ftr2  ■  p0(r,t)  . 


(24) 


Then  the  total  rate  (s  '  )  at  which  particles  exit  the  sphere  r  =  R  is  given  by 


8  r 

Fo(t)  =  ~—  I  dr-4ftr-  ■  p0(r, t)  . 
dt J 


(25) 


Substituting  (22)  into  (24)  and  carrying  out  the  necessary  operations  yields  an  exit  rate  of 

R2 


fo(0=-t 


-R1  IADI 


4ft'/2Di/2t5'2 


(26) 


which,  for  long  times,  decays  as  f  . 


8.4  “FIRST-PASSAGE”  EXIT  RATE  THROUGH  r  =  R 

In  our  Monte  Carlo  model,  when  a  photon  exits  the  sphere  with  radius  r  =  R,  it  is 
removed  from  the  system  and  cannot  return  to  the  inside  of  the  sphere.  Thus  we  model 
the  “first-passage  time”  through  r  =  R  [17]  and  the  corresponding  first-passage  exit  rate. 
The  first-passage  particle  density  inside  the  sphere  is  denoted  p\{r,t).  As  discussed  in 
[17],  this  is  the  solution  to  the  diffusion  equation  (20)  when  an  absorbing  barrier  is  placed 
at  r=R:  the  appropriate  boundary  condition  at  R  can  be  shown  to  be  p\(R,t)=  0.  We  find 
p\(r,t )  using  classical  separation-of  variables.  We  have 

^=DV2p  =  -|^(rp)  .  (27) 

dt  r  or 

We  define 

P\(rd)  =  U  (r)T(t),  U(R)  =  0  .  (28) 
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Then 


1  dT 
DT  dt 


1  d  ,  Tr\  1  2 

- T-trU )  =  const  =  -k 

Ur  dr 2 


The  solution  for  T  is 


rri  -Dk2t 

T  ~e 


(29) 


(30) 


We  set  w(r)  =  rU(r).  Then 


i  /  v  H'  —  v/ 

dr 2 

u  (r )  =  X  (4,  sin(Ar)  +  cos(Ar)) 


(31) 


Since  w  — »  0  as  r  — »  0,  f?„  =  0  for  all  /?.  Since  U{ R)  =  0,  u( R)  =0;  thus  =  nn. 


and 


Pl(rd)  =^  — sin 


A  .  r 


r 


n= 1 


v  f?  y 


exp 


-Dt 


2\ 


(32) 


To  compute  the  4„,  we  note  that 


\dr 


r  ■  sin 


k  nmr  3 


V  f?  J 


A  0%  °)  =  -  Z  4,  J  du  sin(mw)  sin(mi)  =  —  Z  4A,„  ‘  W  =  W  4 

7Z  n  2T  i  2.  2. 


«=1  o 


R 

n  -  i 


n  _  R 


9 


r  ■  sin 


^  m;zr  ^ 

v^y 


A (r, 0)  =  —  •  —  f 4;rr : 2 dr  K  KJs{r )  = 

4nR  R  J  ( mirr\ 


sin 


k  nmr  3 


R 


k  nmr  3 

V  R  y 


m 

2^ 


Then 


,  ,  1  . 
A('%,)=5F§7sm 


k  nnr  3 


v  *  y 


exp 


-At 


v 


r  nn  ^ 

,~R 7 


2\ 


J 


We  now  compute  the  first-passage  time  probability  density  Fi. 

n  R  R 

O 


Fl(t)  = - [  4nr2  dr  ■  px{r  ,t)  =  -4n\  r2dr 

dt J  J 


dPi 

dt 


(33) 


(34) 


(35) 


F i(t)  dt  is  the  probability  of  exiting  the  sphere  in  the  time  interval  (/,  t+dt);  Fj  can 
therefore  be  interpreted  as  the  rate  (s-1)  at  which  particles  exit  the  absorbing  sphere  (“exit 
rate”)  at  r  =  R. 
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dp,  1  ^  n  . 
dt  2  R1  tT  r 


(  n,.2.2^  f 


-Dn 


2  oo  3 


nnr 


V  *  J 


-Dn  7i 


2R 4 


-sin 


6  7 


v  R  J 


exp 


R1 


-Dt 


exp 


-Dt 


f  nn^ 

Cr  j 


r  nn  ^ 

l  ~Rj 


2\ 


(36) 


2 tt3D^  3 

Fi(0=— —L,n  exP 

&  n= 1 


f 

f  nn 

\2 ) 

rR  ,  .  | 

f  nnr 2 

-Dt\ 

ar-r-  sm 

V 

l  ~R. 

’  j 

Jo  1 

l  R  J 

(37) 


The  integral  equals  (/?2/h;r)(-l)n+1;  thus 


r2n  ® 


F.('>=^Z(-l)“Vexp 
A  „=l 


-Dt 


r  nn  ^ 

■  ~R, 


2\ 


(38) 


For  long  times,  the  n  =  1  term  dominates,  and  F/(t)  decays  exponentially  as 
cxpi-D^-2 t/A2).  This  should  be  contrasted  with  the  comparatively  slow  if5'2)  power-law 
exit  rate  (26),  corresponding  to  the  case  where  photons  are  never  removed  from  the 
system. 

The  intensity  (m  “  s  )  of  photons  at  R  is 

7(0  =  -^r^(0  (39) 

4  nR~ 


We  also  note  that  for  small  angle  scattering  the  diffusion  constant  becomes 


D  = 


cL 

3(1- <  cos  6*  >) 


cL  cL  2  cL 


f.  .  e2  "i 

fi  <*2>T 

3  <  e2  > 

3 

1- 

l  2  j 

V  2 

(40) 


and  the  intensity  becomes 


m= 


NncL 

3  ra  <e2> 


^(-!)"+V  exp 


(  nn^ 

l 

(N 

rs _ 

l" ~Rj 

1  3<02  > 

(41) 


This  is  somewhat  different  from  the  corresponding  equation  derived  for  an 
astrophysical  situation  by  Alcock  and  Hatchett  [19];  the  differences  are  discussed  in 
Appendix  B. 


8.5  TELEGRAPH  EQUATION 

The  speed  of  light  does  not  appear  in  the  diffusion  equation  (20),  and  the  solution 
implies  that  a  small  number  of  photons  travel  faster  than  c,  which  is  clearly  non-physical. 
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This  situation  may  be  remedied  by  use  of  the  telegraph  equation  [20],  which  combines 
the  diffusion  and  wave  equations: 


P  = 


1  dp  1  d2 p 
D  dt  c  dt 2 


(42) 


Morse  and  Feshbach  ([16],  p.868)  give  the  Green  function  solution  as 
pOT(r,0)  =  S(r),  pOT(cc,t)  =  0 


- 

8  nD2X 


'  cX' 

v2 Dy 


(43) 


where  I\(x )  is  the  modified  Bessel  function.  (Morse  and  Feshbach  use  a  different 
nonnalization.  To  normalize  so  that  the  integral  over  all  space  is  unity,  their  solution 
must  be  multiplied  by  l/4;zD). 

The  most  complete  theoretical  treatment  of  photon  transport  is  the  Elastic 
Boltzmann  Transport  Equation.  Detailed  solutions  to  this  equation  have  been  obtained  by 
Cai,  Alfano,  et  al.  [21]. 


8.6  SUMMARY  AND  COMPARISON  WITH  MONTE  CARLO  MODEL 

To  summarize:  based  on  the  diffusion  equation  for  isotropic  scattering,  in  a 
uniform  scattering  medium  for  an  impulse  of  photon  emission  at  t  =  0,  r  =  0,  at  a  given 
range  r,  for  long  times  the  density  decays  as  f3/2,  the  total  photon  exit  rate  through  a 
sphere  decays  as  f512,  and  the  first-passage  exit  rate  through  a  sphere  decays  as 
exp (-Dnt/R1)  =  cxp(-/T cLt/OR1))  =  cxp(-/r2 ct/(3R(R/L)).  Figure  18  compares  the 
probability  density  of  first-passage  time  from  diffusion  theory  (solid  lines)  and  from  the 
results  of  the  Monte  Carlo  model  (circles).  The  abscissa  indicates  the  time  after  the  laser 
pulse  strikes  the  target.  Note  changes  in  scales  from  graph  to  graph.  We  again  assume 
that  R  =  100  m.  The  vertical  dashed  line  indicates  the  ballistic  photon  flight  time  =  RJc  = 
0.33  psec. 

The  figure  shows  that  the  results  from  the  two  theories  converge  as  RJL  increases 
and  as  g  decreases,  which  is  expected.  We  conclude  that  the  Monte  Carlo  model  is 
consistent  with  diffusion  theory.  For  the  very  early-arriving  photons  analyzed  in  Section 
6,  the  model  applies  but  the  diffusion  equation  does  not. 
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FPT  Probability  Density  (^is'  ) 


Time  (|is) 

Figure  18.  Comparison  of  Results  from  Diffusion  Theory  and  from  Monte  Carlo  Model 
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9.  CONCLUSION 


The  IDA  Flash  Lidar  Monte  Carlo  model  can  provide  a  useful  tool  for  estimating 
the  quality  of  3D  imagery  from  flash  lidar  through  obscurants  and  for  sorting  out  the 
effects  of  ballistic,  diffuse,  and  “noise”  photons.  In  principle,  moving  targets  could  be 
imaged  if  the  sensor/processor  could  track  the  target,  estimate  its  3D  translation  and 
rotation  as  a  function  of  time,  and  place  photons  in  the  correct  voxels  accordingly.  The 
Geiger  Remapping  Technique  can  be  a  useful  tool  for  interpreting  images  made  using 
Geiger  detectors. 
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APPENDIX  A:  METHOD  OF  GENERATING  RANDOM 
SCATTERING  ANGLES  ACCORDING  TO  THE  SINE-WEIGHTED 
HENYEY-GREENSTEIN  DISTRIBUTION  [7] 


For  Monte  Carlo  simulations  it  is  necessary  to  select  ^-values  from  the  sine- 
weighted  Henyey-Greenstein  (HG)  function,  given  by 

W(6)  =  \  sinflcb(fl)=  ^A-g2)  .  (Al) 

2(l  +  g2  -2gcos6lj 


Here,  the  solid-angle  factor  of  sin  9  is  included  directly  in  the  probability  density  function 
(cf.  Eq.  2)  and  the  normalization  condition  is 


In  n  n 

J  d(j)\w(9)d9  =  2n\ 


sin#(l -g)d6 
2(l  +  g2  -2gcos(?) 


=  1 


(A2) 


(The  g  =  0  case,  namely  isotropic  scattering,  is  easily  checked.) 


We  now  wish  to  generate  random  values  of  6  with  probability  density  given  by 
(Al).  A  general  method  for  generating  random  numbers  with  a  specified  density  is 
outlined  in  [8],  Chapter  7.  Let  x  be  unifonn  in  [0,  1];  that  is,  x  has  density  function 


p{x)  = 


0  <  x  <  1 

otherwise 


(A3) 


Assume  there  is  a  function  0=  /(x)  that  transfonns  p  into  W.  By  conservation  of 
probability,  W  must  satisfy 

\p(x)dx\  =  \W(0)dO\  .  (A4) 


Rearranging,  one  finds 

Te~~w ^ 


which  is  easily  integrated  to  yield 

x=\w(0)d6=F(O ) 


(A5) 


(A6) 
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Here  F(9)  is  the  cumulative  distribution.  The  desired  transformation  is  given  by  the 
inverse  of  the  cumulative  distribution  function,  namely 

d=f{x)  =  F-\x )  (A7) 


In  order  to  arrive  at  a  closed-form  solution  for  f  we  require  both  a  closed-form 
expression  for  the  cumulative  distribution,  and  a  closed-form  expression  for  the  inverse. 
For  the  sine -weighted  HG  function,  both  closed-form  expressions  exist. 


Using  equation-solving  software,  the  cumulative  probability  function  is  found  by 
computing  the  integral 


F<m=\ 


r  sin£(l-g2)d£  _(l-g2) 

r  i 

1 

o  2^1  +  g2  — 2gcos4r)  ^g 

[i -g  , 

Jl  +  g2 -2gcos0  y 

(A8) 


and  the  inverse  function  is  found  to  be 


0=f(x)  =  F  1(x)  =  cos 


l-2x  +  2g(l  -  x  +  x^ )  +  g2  (1  -  2x)  +  2g\x(x  - 1) 
(2gx-l-g)2 


.  (A9) 
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APPENDIX  B:  MODEL  OF  ALCOCK  AND  HATCHETT— SMALL- 

ANGLE  SCATTERING 


Using  the  equations  of  radiative  transport  in  a  scattering  medium  and  under  the 
assumption  of  very  small  scattering  angles,  Alcock  and  Hatchett  [19]  computed  an  exact 
solution  for  photon  scattering  vs.  angle  and  time.  Their  results  have  been  frequently 
quoted  and  used,  for  example,  in  [22].  A  brief  summary  of  their  results  follows. 

In  the  notation  of  [19],  <j>=  angle  between  photon  trajectory  and  radius  vector,  0  = 
scattering  angle,  k=  scattering  coefficient  =  ML,  r=  optical  depth  =  R/L,  N=  total  photon 

number,  y  =  ct/L,  and  £=  (f> .  Alcock  and  Hatchett  obtain  an  expression  for  the  intensity 

—2  —1 

(m  s  ),  given  by 

A J 

I(M  =  -T-IP(g,r,T)  .  (Bl) 

4  n-R- 


In  the  limit  of  t  — »  oo,  that  is,  a  large  number  of  scatterings  they  show  that 


P(C,y;r)  = 


1 


<  6-  >2  r3 


(T 


Y 


T  <6  >  T  <0  > 


(B2) 


where  G  is  a  complicated  integral  function  given,  and  graphed,  in  [19].  The  integral  over 
angle  is 


P(y;r) 


4ft2 

t2  <  e2  > 


]^(-l)"+1«2exp 


n= 1 


v 


-In2  ft2 y  ^ 
t2  <  6 2  >y 


(B3) 


—2  —  1 

In  our  notation,  the  intensity  (m  s  )  vs.  time  is  found  from  (Bl)  and  (B3)  to  be 
,, ,  NcL 


R*  <  e2  >  f- 


X(-ir‘»!exp 


=1 


f  nft 'j 

i2  2 cLt 

A 

to 

V 

1 _ 

(B4) 


By  comparison,  our  (43)  is 


m= 


NftcL 
3 R4  <02> 


X(-l)"*'n!exp 

n= 1 


f  nft'}, 

l 

(N 

<N 

l  ~R) 

1  3<02  > 

(B5) 


The  coefficient  of  (B5)  is  ft/3  times  that  of  (B4);  more  importantly,  the  exponent  of  (B5) 
is  1/3  that  of  (B4).  These  differences  are  discussed  in  Appendix  C. 
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APPENDIX  C:  DIFFUSION  CONSTANT 


This  Appendix  documents  the  construction  of  appropriate  rotation  matrices  with 
which  to  implement  the  random  walk  for  the  Ballistic  Imaging  model.  In  the  limit  of  a 
large  number  of  scattering  events  the  random  walk  can  be  treated  as  an  isotropic  diffusion 
problem;  the  diffusion  constant  for  this  limit  is  derived 

By  convention,  the  photon  travel  direction  is  chosen  to  be  along  the  positive  z  axis 
of  the  current  frame  (the  “photon  frame”).  That  is,  in  the  current  photon  frame  the 
direction  of  travel  is  always  (0,0,1).  The  net  photon  displacement  must  be  computed  in  a 
global  coordinate  system  (“lab  frame”).  The  transformation  from  photon  frame  to  lab 
frame  is  accomplished  by  a  suitable  product  of  rotation  matrices,  developed  below. 

Each  scattering  event  is  described  by  two  angles,  9  and  (/>.  The  polar  angle  6  is 
measured  from  the  positive  z  axis  of  the  frame  before  scattering  to  the  positive  z  axis  of 
the  frame  after  scattering,  and  is  chosen  from  the  sine -weighted  Henyey  Greenstein 
probability  density.  The  azimuthal  angle  (j)  is  measured  from  the  positive  x  axis  of  the 
frame  before  scattering  to  the  positive  x  axis  of  the  frame  after  scattering,  and  is  chosen 
from  a  uniform  probability  density  on  [0,  2n).  The  path  length  between  two  scattering 
events  is  selected  from  an  exponential  probability  density 

f{l)  =  yLe,,L  ,  (Cl) 


where  L  is  the  mean  free  path. 

Figure  C-l  shows  the  relation  between  the  frames  before  and  after  a  scattering 
event.  The  (x,  y,  z)  frame  is  before  scattering,  while  the  {x",y",z")  frame  is  after  scattering. 
The  rotation  order  convention  is  to  rotate  first  by  (j)  about  the  z  axis  to  yield  the  (x\  y' ,  z') 
frame,  then  by  #  about  the  x'  axis  to  give  the  (x",  y",  z ”)  frame.  The  net  rotation  is  built  up 
in  two  steps.  The  coordinates  of  a  vector  in  the  (x',  y\  z')  frame  are  expressed  in  terms  of 
the  coordinates  in  the  (x,  y,  z)  frame  through 

^cos  (f)  -sin^  0N 

sin^  cos  (j)  0  ]|  y'  \  =  Rz(<f>)\  y'  \  .  (C2) 
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Figure  C-1.  Coordinate  Systems 

Similarly,  the  coordinates  of  a  vector  in  the  (pc',  y',  z')  frame  are  expressed  in  terms  of  the 
coordinates  in  the  (x",  y",  z")  frame  through 
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The  net  rotation  relating  coordinates  before  and  after  scattering  is 
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The  nth  scattering  event  is  described  by  the  rotation  matrix  Mn=M(<f>n,0n).  Since  the  photon 
direction  is  always  (0,0,1)  in  the  current  photon  frame,  the  photon  direction  in  the  lab 
frame  is  found  by  transforming  (0,0,1)  backwards  n  times  via  the  matrix  product 


vn  =  MXM2...M  h 


0 


(C5) 


Letting  Sn  denote  the  matrix  product  at  the  nth  scattering  event,  and  k  denote  (0,0,1)  in 
the  current  photon  frame,  one  then  arrives  at  the  “scattering  direction  map”  given  by 

v  =51  k 


Sn+,  =  SnMn+l  ;  50  =  1 


(C6) 


In  the  lab  frame  the  net  displacement  of  the  photon  after  n  steps  is 
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(C7) 


n  n  n 

R„  =  Zr «  =  Z  L^a  =  Z(Ak 

a=  1  a=  1  a=  1 

where  la  is  the  step  length  at  step  a. 

The  diffusion  constant  for  the  random  walk  is  found  as  follows.  The  mean  square 
displacement  after  n  steps  is 

v 

kj 


where  T  denotes  the  transpose  operation  and  <>  denotes  ensemble  average.  Simple 
manipulation  yields 

<K1>=t</«1kTS„Xk>  +2£  £  <‘JfkT  SjSfk>  .  (C9) 

a=\  a= 1  (3=a+\ 

The  first  sum  (diagonal  terms)  evaluates  to  n  </2>,  since  the  Sa  are  orthogonal 
matrices  and  satisfy  Sa  =  Sa  .  In  the  second  sum  (off  diagonal  terms)  the  average  over 
step  lengths  and  angles  factors,  since  the  step  lengths  and  the  scattering  angles  come  from 
statistically  independent  physical  processes.  Also,  the  steps  lengths  are  independent,  so 
</a  /p>  =</„></p>  =  </>2.  Applying  these  facts  to  (C9)  yields 

<Rn2  >=n</2>  +2</>2  Zkr<5/^>k  .  (CIO) 

/3>a 


The  average  in  the  sum  is  given  by 

<  SaTSp  >  =  <  (M  M2  ...Ma)T(Mx  M2  ...Mp)  > 
=  <  Mj Ma_y  ...  MlTMlM2...Mp  > 

=  <K+lMa+2-Mp> 

=<M>p~a 


(Cll) 


where  the  third  equality  follows  from  orthogonality  and  the  last  equality  follows  from  the 
independence  of  the  scattering  events.  Using  the  fact  that  <siru/J>  =  <cos0>  =  0  the 
average  rotation  matrix  is  found  from  (C4)  to  be 

"0  0  0  ^ 

<  M  >=  0  0  0  .  (Cl 2) 

v0  <sin<9>  <cos6>y 

Therefore 
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0 


and  it  follows  that 


0  0  ^ 
<  M  >k  =  0  0  0 

0  <  cos<9  >A '<  sin<9  >  <  cos 0  >k 

<R2>  =  n<l2>  +  2  </>2  EE  <  cos  6  >p  “ 

a= 1  fi=a+l 


(Cl  3) 


(Cl  4) 


The  double  sum  can  be  evaluated  by  recognizing  there  are  (n-k)  tenns  of  <cos8>k,  hence 

]T  £  <  cos 9  >p~a  =  (w  -  fc)  <  cos 6  >k  .  (C15) 


a= 1  p=a+\ 


k= 1 


The  second  sum  in  (Cl 5)  can  be  evaluated  by  standard  techniques1  and  the  result  is 

£4.  ,  k  <cos8  >n+l -n  <co$0>2 +{n  -\)<  co$0> 

Y^in-k)  <cosd>k  = - - - — - - -  .  (Cl 6) 

(l-<cos#>) 

Since  the  steps  are  chosen  from  an  exponential  probability  density  the  mean  square  step 
length  and  mean  free  path  are  related  via2  <I2  >  =2<l  >2=2L2 .  The  mean  square 
displacement  is  then 


<R2  >  =  </2  > 


n  + 


<  cos  0>{<  cos  8  >"  -n  <  cos  8>+n- 1 ) 
(1-  <  cos#  >)2 


(Cl  7) 


For  large  n  and  for  sufficiently  large  scattering  angles  6  (to  be  discussed  further  below), 
<  cost?  >"  can  be  neglected  in  the  numerator  of  (Cl  7)  and  one  arrives  at 


<  R,  > 


n<l2  > 
l-<  cos  6*  > 


(CIS) 


For  a  three-dimensional  random  walk  the  diffusion  constant  is  related  to  the  mean  square 
displacement  and  the  flight  time  through  [23] 


n- 1  i  n-l 

First  rewrite  the  sum  as  <  COS  8  >"  '  =  —  T"!  id' 


1  d 


f  H- L  A 


J=  1 


7  7=1 


n— 1  o 

q  oq 


J  ,  where  q  =  (<cos  6>) 


VM  J 


The  sum  inside  parenthesis  is  a  geometric  sum  that  evaluates  to 


q-q 
i —q 


.  Carrying  out  the  necessary 


operations  yields  the  expression  in  (Cl 6). 

00  -.00 

The  mean  square  step  length  is  given  by  <  /  >=  /  p(l)dl  = —  /"  e  /L  dl  =2 L~  . 

L 
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<R 2  >=6Dt 


(Cl  9) 


Combining  (Cl 8)  and  (Cl 9)  and  using  the  fact  that  after  n  scatterings  the  average  path 
length  is  nL=ct,  where  c  is  the  speed  of  light  and  that  <  I2  >  =  2Lr  ,  the  diffusion 
constant  is  given  by 

D  = - — -  .  (C20) 

3(1-  <  cos  6>) 

We  now  examine  the  restriction  on  the  scattering  angle  6  in  arriving  at  the  approximation 
in  (Cl 8).  The  approximation  starts  with 

<  cos  6  >"«  nl(-  <  cos  9  >)  ,  (C2 1 ) 

which  for  small  angles  becomes 

<  1  -£->"«  n{\-  <  4-  >) 

1-f  <92  >«f  <92  >  ,  (C22) 

n  <6~  >  »  1 

or 

M-m,  »  1  •  (C23) 

The  central  assumption  in  Alcock  and  Hatchett’s  treatment  [19]  of  X-ray 
scattering  is  that  of  extremely  small  scattering  angles  and  a  relatively  small  number  of 
scattering  events,  characterized  by  the  opposite  condition  to  (C23),  namely  [24] 

M„,  «  1  •  (C24) 

The  diffusion  analysis  outlined  in  this  Appendix  does  not,  therefore,  apply  to  their 
treatment.  In  fact,  Alcock  and  Hatchett  say  (pp.  459-460),  “These  formulae  are  all 
developed  under  the  assumption  that  the  resultant  angle  of  the  photon  trajectory  after  any 
number  of  scatterings  is  small.  This  requirement  keeps  the  analysis  in  a  quantitatively 
different  regime  from  that  of  random  walk  analyses,  and  the  results  are  qualitatively 
different.” 
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